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Abstract The evolutionary effect of recombination depends crucially on the epistatic in- 
teractions between linked loci. A paradigmatic case where recombination is known to be 
strongly disadvantageous is a two-locus fitness landscape displaying reciprocal sign epista- 
sis with two fitness peaks of unequal height. Although this type of model has been studied 
since the 1960's, a full analytic understanding of the stationary states of mutation-selection 
balance was not achieved so far. Focusing on the bistability arising due to the recombina- 
tion, we consider here the deterministic, haploid two-locus model with reversible mutations, 
selection and recombination. We find analytic formulae for the critical recombination prob- 
ability Vc above which two stable stationary solutions appear which are localized on each 
of the two fitness peaks. We also derive the stationary genotype frequencies in various pa- 
rameter regimes. In particular, when the recombination rate is close to and the fitness 
difference between the two peaks is small, we obtain a compact description in terms of a 
cubic polynomial which is analogous to the Landau theory of physical phase transitions. 

Keywords evolution of recombination • reciprocal sign epistasis ■ bistability ■ Landau 
theory 
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1 Introduction 

After more than a century of research, the evolutionary basis of sex and recombination re- 
mains enigmatic lOttoL l2009h . In view of the complex evolutionary conditions faced by 
natural populations, the search for a single answer to the question of why sexual reproduc- 
tion has evolved and is maintained in the vast majority of eukaryotic species may well be 
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Fig. 1 Schematic of a two-allele, two-locus fitness landscape with reciprocal sign epistasis and unequal peak 
heights. Fitness is represented by the height above the plane in which the four genotypes (labeled 0, 1, 2, 3) 
reside. 



futile. Nevertheless, theoretical population genetics has identified several simple, paradig- 
matic scenarios in which the conditions for an evolutionary advantage of sex can be identi- 
fied in quantitative terms, and which are therefore open (at least in principle) to experimental 
verification. 

One such scenario was proposed in the context of the adaptation of a population in a 
constant environment encoded by a fitness landscape, which assigns fitness values to all 
possible genotypes. In this setting, the relative advantage of sexual reproduction with re- 
spect to, say, the speed of adaptation or the mean fitness level at mutati on- selection balance, 
depends crucially on the e pistatic interactions between different loci dde Visser and ElenS . 
I2OO7I : IKouvos et all l2007h . In its simplest form, epistasis is associated with the curvature 
of the fitness surface, that is, the deviation of the fitness effect of multiple mutations, all of 
which are either deleterious or beneficial, compared to that predicted under the assumption 
of independent mutations which contribute multiplicatively or additively to fitness. It is well 
established that recombination speeds up the adaptation in such a fitness landscape if the 
epistatic curvature is negative, in the sense that the fitne ss of multiple mutants is lower than 
expected for independent mutations dKondrashovl . 1 1 98 sh . 

Unfortunately, experimental evidence indicates that negative epistasis is n ot sufficiently 
widespread to provide a general explanation for the evo lution of recombination dde Visser et all 
1 19971 : Elena and Lens kil, ll997l : lBonhoeffer et a]ll2004[) . Instead, empirical studies have high- 
lighted the importance of a more complex form of epistasis, termed sign epistasis, in which 
not only the magnitude but also the sign of the fitness effect of a given mutation depends 
on the presence of mutations at other loci jWeinreich et alll2005l ). A simple example of sign 
epistasis is shown in Fig.[U which depicts a haploid two-locus fitness landscape. In this case 
the sign epistasis between the two loci is reciprocal, wh i ch lea ds to the appearance of two 
fitness peaks separated by a fitness valley jPoelwiik eta]Ll2007h . 

In the two-locus setting sign epistasis can be viewed as an extreme case of positive 
epistasis, and it may be expected to lead to a disadvantage of recombination. In a recent 
simulation study of the effect of recombination on empirically motivated five-locus fitness 
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landscapes displaying sign epistasis jde Visser et a il2ooi), we found that recombination 
can hamper adaptation on such a landscape in a rather dramatic way: Instead of moving to 
the global fitness optimum, populations get trapped at local optima from which (in the limit 
of infinite population size) they never escape. Mathematically, these numerical calculations 
suggest the appearance of multiple stable stationary solutions of the deterministic, infinite 
population dynamics above a threshold value of the recombination rate. 

The simplest situation where this phenomenon can occur is the haploid two-locus land- 
scape shown in Fig. [TJ where it implies a bistability with two stationary solutions local- 
ized near each of the fitness peaks labeled by and 3. Indications for the occurrence of 
bista bility in this system can be found in several earlier studies of the two-locus prob- 
lem. |c row and Kimurj (Il965h derived a condition for the initi al increase of the high peak 
mutant in a population dominated by the low peak genotype. lEshel and Feldma J( ll970f) 
established a sufficient condition for the low peak population to remain trapped for all 
times when mutations are unidirectional (0 -^h 2 — >■ 3), which w as sub sequently refined 
bv iKarlin and McGreeod ( Il97lh . Ipeldmanl ( Il97lh and Ikutschmarj Jl994h obtained condi- 
tions for the existence of multiple stationary solutions in the absence of mutations. The case 
of symmetric fitness pe aks was conside red a s a model for compensatory mutations in RNA 
secondary structure bv IStephar? (1996) and Higgs' (^ 1998h . who also addressed the role of 
genetic drift. Recent studies have considered the effect of recombination on the dynamics of 
peak escape in the asymmetric two-lo cus landscape for finite as well as infinite populations 
( IWeinreich and Chao[|2005l : |jahH 2010). 

The corresponding diploid problem has also bee n studied extensively ( lKimural,[l95^ : 
iBodmer and FelsensteinL [l967h . iBtiree] ( Il989l . [2OO0I) established conditions for bistability 
in a class of quantitative genetics models of stabilizing selection that are formally equiva- 
lent to the diploid version of our problem with symmetric fitness peaks. Finally, bistability 
induced by rec ombination has been observed in multilocus mo dels in the context of quasis- 
pecies theory teoerl ijst et aj 1 19961 : IJacobi and NordahlL l2006h and evolutionary computa- 
tion IWright et al, 2C^. 

However, a comprehensive analysis of the paradigmatic haploid two-locus system with 
reversible mutations, reciprocal sign epistasis and fitness peaks of unequal height does not 
seem to be available, and the present paper aims to fill this gap. In the next section we intro- 
duce the evolutionary dynamics used in this work. We then show that finding the stationary 
solutions of the model amounts to analyzing the zeros of a fourth order polynomial, and de- 
vote the bulk of the paper to extracting useful information about the critical recombination 
rate, the genotype frequency distribution, and the mean fitness in various parameter regimes. 
We conclude with a summary of our results and a discussion of open problems. 



2 Model 

We consider a haploid, diallelic two-locus system. The alleles at a locus are denoted by and 
1, and the four genotypes are labeled according toOO— >-0, 01— >-l, 10— 5'2, 11— >-3. Each 
genotype ; is endowed with a fitness w,, which defines the fitness landscape. The population 
evolves in discrete, non-overlapping generations under the influence of selection, mutation, 
and recombination. Since we are interested in the emergence of bistability, in the sense of the 
existence of multiple stationary frequency distributions, the limit of infinite population size 
is assumed. Thus, the frequency of each genotype evolves deterministically. The frequency 
changes according to the following order: selection, mutation, then recombination. 
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Let fi{l) denote the frequency of genotype / at generation T. The frequency of the geno- 
type at generation T + 1 after selection is proportional to the product of its frequency at 
generation T and its fitness w/. After selection, mutations can change the frequency of geno- 
types. We assume that the mutation probability per generation and locus, /i, depends neither 
on the location of the locus in the sequence nor on the alleles, and mutations are assumed 
to be independent of each other. Mathematically speaking, the mutation from sequence j to 
sequence / then occurs with probability 

[/,.^ = (i-^)2-"'ai)/M^ (1) 

where d{i, j) is the Hamming distance between two sequences / and j, i.e. the number of loci 
at which the two sequences differ. After selection and mutation, the frequency distribution 
will be 

l{r+i)=l^U<j^^fj{T), (2) 

where w(t) = T.i^ifii'^) is the mean fitness at generation t. 

After selection and mutation, recombination follows. At first, two sequences, say j and 
k, are selected with probability fjfk- With probability 1 — r, no recombination happens. With 
probability r, however, the two sequences recombine in a suitable way and generate a recom- 
binant, which may be identical to or different from j and k depending on the recombination 
scheme and the properties of the two initial sequences. One may interpret the case r < 1 as 
a model for organisms which can reproduce sexually as well as asexually. We will use the 
one-point crossover scheme, which is summarized as 

with prob. r/2, 
with prob. r/2, 

with prob. (1 -r)/2, ^ 
with prob. {I — r) /2, 

where Xi (X2) is the allelic type at the P' (2"'') locus, j = X1X2, k = x'yXi ^^e parental 
genotypes, and the four types on the right hand side are the possible resulting recombinants 
with their respective probabilities. If denotes the probability that the resulting allelic 
type is i in case types j and k (not necessarily different) attempt to recombine, the frequency 
of type i at generation T-|- 1 becomes /J = T,jkRi\jkfjfk- 

After selection, mutation, and recombination, the frequency of each sequence in the next 
generation becomes 
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where 5(t) 


= /0(T)/3C 


r)wo>V3 — /i (t)/2(t)m'iW2 and /i's and fl's mean the frequency at 



considered in the literature corresponds to r = ^ in Eq. Q. 
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By rearranging (|4]i, we get 



^^(/0-/3) =/0W0-/3W3, (5) 

(fi ~ fi) = /i - /2W2, (6) 

P(T) 



A nontrivial conclusion one can draw from Eq. ^ is that linkage equilibrium, /q/j = /1/2, 
is attained in one generation if r = 1, which corresponds to the one-point crossover scheme 
with recombination probability 1 . 

The stationary distribution is calculated by setting /J = // in the above three equations, 
which gives 



(8) 
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B, (10) 
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where w is the mean fitness at stationarity. With the two additional conditions 

/0 + /l+/2+/3 = l (11) 

and 

^ = ^0/0 + 1^1/1-1-1^2/2 + 1^3/3 (12) 

the steady state solution is fully determined. Without loss of generality, W3 is set to be 
largest (with possible degeneracy). For simplicity, we set wi = W2, which by ([9]l implies that 
/i = /2 = /. Since the dynamics is invariant under multiplication of all fitnesses by the same 
factor, we can choose 

1^3 = 1, Wo = 1 — f, wi = 1V2 = 1 — f — .? (13) 

with < f < 1 and —t<s< \ —t. For the sake of brevity, however, we will sometimes 
use Wo and wi rather than s and t in what follows. The behavior for unequal valley fitnesses 
(wi 7^ W2) will be discussed in Sect l3.8l 

Note that the problem studied by iHiggsl ( Il998h corresponds to the case with t = and 
< i. In this paper, we will assume that t > Q and j > 0, that is, the fitness landscape has 
a unique global optimum and reciprocal sign epistasis (Fig.O. Unlike the fitness landscape 
with symmetric peaks, it is difficult to find the solution exactly (see also Appendix iBt. We 
will investigate the approximate solutions by assuming that some of the parameters are very 
small compared to others. 
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3 Bistability 

3.1 General behavior of solutions 

Since the frequency of each genotype is strictly positive in the steady state, Eq. ^ excludes 
the mean fitness w from being in the range between (1 — 2n)wo and (1 — 2/i)w3. For obvi- 
ous reasons, we will refer to a solution with w > (1 — 2/i)w3 as high-fitness solution (HFS) 
and a solution with m> < (1 — 2jU)wo as low-fitness solution (LPS). As an immediate conse- 
quence of Eq. we see that a HFS (LPS) implies /; > /o (/? < /o). The largest fitness 
being W3, the existence of a HPS is expected regardless of the strength of the recombination 
probability (we will see later that this is indeed the case). Hence we are mainly interested in 
the conditions which allow for a LPS. Later, we will see that if it exists, there are actually 
two LPS's, only one of which is locally stable. Since we are interested in stable solutions, 
in what follows LPS refers exclusively to the locally stable solution. Intuitively, the HPS 
should be locally stable, so the emergence of a LPS naturally implies bistability. 

Without much effort, one can easily find necessary conditions for the bistability. Pirst 
note that w is larger than wi by defirution. Therefore, a LPS is possible only if wi < (1 — 
2/i)wo, or 

^<2(T^- (14) 

If we put jU = I in Eq. (|4]l, every fj becomes ^ regardless of r and the /'s. Since this is a 
unique equilibrium state for /x = 5 , it could have been anticipated that bistability requires a 
restriction on fi. A necessary condition on r can be obtained from Eq. dlOb . While the nu- 
merator of the expression defining B is always positive, the denominator would be negative 
for the LPS if wq — (1 — r)w3 < 0. Hence a necessary condition for bistability is 

r>l-^=t, (15) 
W3 

which appears also in earlier works jCrow and KimuraL 1 19651 : lEshel and FeldmanL Il970l : 
iKarlin and McGregoil[T97ll : IPeldmarl 1 97 iL iRutschman. 1 1 994) . Most of the calculations in 



this paper are devoted to refining the conditions for bistability. To this end, we will reduce 
the five equations ( 1819110111112b to a single equation for w. 

Equations ^ and dlOb along with the normalization jl lb yield the relations 

^ = 2v^ + V^(l+A)' ^^ = 2v^+v^(l+A)' (''^ 

where we have used the fact that the /, 's should be positive, and the definition ( 112b of the 
mean fitness w implies that 

2{w-wi) = \ -{wi + WQA-{\+A)w). (17) 
V ^ 

Since our main focus is on the LPS, it is convenient to define the auxiliary variable x through 

w=(1-2m)(wo-x), (18) 

which implies that x < —t and x > for the HPS and the LPS, respectively. Note that x is 
equivalent to the mean fitness and equilibria will be found in terms of x. We also note for 
future reference that with this reparametrization, A and B in Eqs. ^ and dlQb become 

A=l + ; B=X + {l-r)- ;r.7^, (19) 

X (Wq — xy — {\ — rjWQWi 
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Taking the square on both sides of Eq. ( 117b results in a quartic equation 

h{x) = ho{x) + rhi [x) = 0, (20) 

where the polynomials hQ{x) and h\ {x) are independent ofr (see Appendix lAl and the Math- 
ematica file in online supplement for explicit expressions). 

We can draw some general conclusions concerning solutions of the quartic equation by 
evaluating h{x) at selected points. Since h{x) is negative at x = 0, x = —t, and at the point x = 
xi defined by (1 — 2pL){wQ — x\) = w\, and the coefficient of fourth order term is positive (see 
Appendix|All, there always exist solutions in x > x\ and x<—t which correspond to vt> < wi 
and vi> > (1 — 2/i)w3, respectively. The meaningless solution, w <w\, has appeared because 
we took squares on both sides of Eq. | |17| |. In the Mathematica file in online supplement, 
we show that h(x) is positive when x= 1— f— 1/(1 — 2/i),orvi> = W3 = 1. This proves, as 
anticipated, that the HPS with the mean fitness in the range (1 — 2jU)w3 < vP < 1^3 is present 
for all values of r. 

Hence the condition for bistability is recast as the condition for h{x] to have a positive 
solution which is smaller than x\. Because /!(0) and h.(x\) are negative, the existence of a 
solution in the range < x < x\, <x equivalently w\ < w < (1 — 2/i)wo, always entails two 
solutions in the same region if we count the number of degenerate solutions (that is, two 
identical solutions) as 2. Let us assume that h{x) = has two degenerate solutions at jc = x^. 
when r = Tc (as we will see, is the critical recombination rate above which bistability 
exists). This means that Xc is the simultaneous solution of two equations h(xc) = h'{xc) = 
0, where the prime denotes the derivative with respect to the argument. Later, this simple 
relation will play a crucial role in finding as well as Xc- Now let us change r infinitesimally 
from rc to + e,., and let the solutions of h{x) =0 for r = r^. + take the form Xc + Ex- Note 
that we are only interested in solutions with e, — )• as e,- — )• in the complex x plane. Since 
two other solutions exist outside of the range < j: < and both h{0) and h{x\ ) are negative, 
h{x) with r = has local maximum at jc = x^, that is, h"{xc;rc) < 0. Figure |2]illustrates this 
situation. 

Using Eq. i l20t we get 

h{xc + £x;rc + er) w -^\h"{xc;rc)\£^ + £rhi(xc)=0. (21) 

Hence real solutions are possible if the second term is positive. By definition, rchi{xc) = 
—ho{xc). For r = 0, Eq. ( 120b reduces to ho{x) = 0, and we may conclude from the condition 
Eq. ( 115b . which is violated for r = 0, that this equation does not have a solution in the region 
< X < xi. Hence /lo(xc) < because ho(p) < and hQ{xi) < (see Appendix lAb. which, 
in turn, implies that hi (xc) > 0. To summarize, we have proved that if Xc is the degenerate 
solution of h{x) =0 when r = r^, there are two solutions in the region < x < xi when 
r > Tc. One should note that r^, if it exists, is unique, as otherwise a contradiction to Eq. ( 121b 
would arise. 

To establish the existence of bistability for r > rc, it remains to find r^. Even though 
general solutions of quartic equations are available, it is difficult to extract useful information 
from them. Rather, we will look for approximate solutions by assuming that one of the 
three parameters /i, s, and t is much smaller than the other two. In fact, it follows from the 
condition ( 114b that there cannot be any bistability when .? <C /x (unless f — >■ 1). Hence, in this 
paper, we will not pursue the case where s is the smallest parameter. 

Before turning to the derivation of the approximate solutions, we further exploit the lin- 
ear r-dependence of h{x} in Eg. (120b. It implies that the two equations h{xc;rc) = h'{xc;rc) = 
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-0.04 ^ ' ' ' ' ^ 

-0.4 -0.2 0.2 0.4 0.6 

X 



Fig. 2 Behavior of the function h{.x) around the critical recombination probability. In this example, i = 0.5, 
; = 0.4, and ^ = 0.01 are used, which yields ~ 0.965. The curves meet at (three) points where h[{x) 
vanishes. The locations of the points x = —t and x = x\ are indicated by filled circles. Low fitness solutions 
correspond to zero crossings in < A' < a'i, and high fitness solutions to those with x < — f. The HFS's are 
indicated by an arrow, which happens to be close to — f. The zero crossings with x > x\ are spurious. Inset: 
Close-up view of the boxed area. For clarity, the vertical line at .v = is also drawn. One of the solutions of 
h\(x) = happens to be close to x = 0. Note that for r = r^, two solutions become identical (degenerate). For 
r> Vc, the LFS is indicated by an aiTow. 



with two unknowns can be reduced to a single equation for Xc which does not involve r^. 
To be specific, from ho{xc) + rchi (xc) = /lo(jCc) + ^ch'^ (xc) = we obtain 

_ hojx,) _ h[^{x,) ^^^^ 
hi{xc) h\{xcY 

or 

H{x,) = h,{x,)h\ {x,) - hi {x,%(x,} = 0. (23) 

Thus, rather than finding and Xc simultaneously, we will first find Xc by solving Eq. ( 123b , 
which in turn, gives from Eq. ( 122b . Equation ( 123b will be analyzed below, after we have 
introduced one more useful concept. 



3.2 Critical mutation probability 

As evidenced by Eq. ( 114b , a finite critical recombination probability can exist only if jU 
is sufficiently small. Mathematically, this implies that r^. diverges as /i approaches a certain 
critical mutation probability fie, such that bistability is not possible for /i > /ic- Although r 
cannot, strictly speaking, exceed unity, we will see later that this definition of will be of 
great use to find an accurate expression for Kc- Setting = °° in Eq. ( 122b . we see that is 



9 



the solution of the equations hi (xc) = h\ {xc) = 0. Since hi {x) is a cubic function taking the 
form hi {x) = —Ct,x^ — C2X^ +Cix — Cq, Xc is also the solution of the equation 

G[x) = xh[ (x) - 3hi (jc) = C2X^ -2Cix + 3Co = 0. (24) 

From G{x) and h\ (x), we can construct two linear equations for x such that 

Gi{x)=C2h'i{x)+3C3G{x) = -2(3CiC3 + C|)x + CiC2 + 9CoC3 =0, (25) 

GzW = ^{CiG{x)-3CQh'i{x)) = (CiC2 + 9CoC3)x-2(C?-3CoC2) =0, (26) 

where we have used the fact that Xc ^ 0. Hence, the value of Xc for r,; = °° is given by 

^„ ^ C1C2 + 9C0C3 ^ 2(Cf-3CoC2) 
2(C| + 3CiC3) CiC2+9CoC3' 

Note that C,'s depend on /I, and t, which means we have an equation for jic such that 

(CiC2 + 9CoC3)2-4(cf-3CoC2)(C| + 3CiC3) =0, (28) 

which, in turn, provides x'^ by inserting /ie in Eq. ( 1271 ). In fact, Eq. ( 128b is equivalent to the 
discriminant of cubic polynomials (see the Mathematica file in online supplement). 

We now assume that s,t <C 1, which also implies 1 by Eq. ( I14l l. Then to leading 

order the C, 's become 

Ci=2s + t, C2 = (^ + 2/l)(2J + ?)-^^ Ci = t{s^ -2n{2s + t)), Co = M^, (29) 

and Eq. ( 128 b becomes (see the Mathematica file in online supplement) 

3i"'v2(l - f) (1 + V - z(2 + v)fM{z) = 0, (30) 

where z = pi/s, v = t /s, and M(z) is a cubic polynomial with 

M(z) = 32{v + l)z^ - (13v2 +48v +48)z2 + 2(2v3 +7v2 +9v + 6)z- (1 + v)2. (31) 

Note that z = (1 + v)/(2 + v) cannot be a solution because of Eq. ( 114b . Hence the critical 
mutation probability is the solution of the equation M{z) — 0. As shown in the Mathematica 
file in online supplement, M(z) is an increasing function of z, which, along with M(0) < 0, 
allows only one positive real solution of M(z) = unless v = 0. One can find the exact solu- 
tion in the Mathematica file in online supplement, which is not suitable to present here. We 
will just present the asymptotic behavior of the solution for later purposes (see the Mathe- 
matica file in online supplement). When r <Cs' (v <C 1), 



s 



2/3' 

1-3 ■ ■ ■ 



(i)' 



1/3 

andx~=(^) , (32) 



and when i < f (v S> 1) 



.2 



... 4,and.r = ^. (33) 

Although we derived the above two expressions from the exact solution (see the Mathe- 
matica file in online supplement), it is not difficult to find the asymptotic behavior without 
resorting to the exact solution. When v is finite, it would be more useful to have a numerical 
value. In the case v = 1 (f = i) we get 

\ic w 0. 107j and w 0.0616j. (34) 

Below we will see how yLc can be used to derive improved approximations for Vc- 
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3.3 Approximation for small mutation probability 

Now we will move on to finding the critical recombination probability. We begin with the 
investigation of the approximate solutions for small mutation probability (jU <C s,t). Let us 
assume that = jcq + a^/i + 0(/J^), which should be justified self-consistently. For later 
purposes, we introduce the parameters 

a = {l-t){s + tf-s^, p = {l-t){s + tf + s^. (35) 



Note that a is positive if ^ < —t + 1 — f , which is automatically satisfied because s is 
smaller than 1 — r by definition. 

The leading behavior of H(xc) becomes 

H{xc) = xl{t +xo)^{a2xl + 2aixi) + aQ) + 0(^), (36) 

where the a, 's are parameters satisfying — 00^2 < (see the Mathematica file in online 
supplement). Hence, the solutions for xq are 0, —t, and two complex numbers. Since Xc 
must be positive, the only possible solution for xq is xq = 0. Accordingly, the actual leading 
behavior of H(xc) becomes (contributions of order 0(1) and 0(/i) vanish) 



H{xc) = -li^s^f- ((1 - tfa - alP) + 0(jLi ) = 0, (37) 



«M = (l-0(^) ■ (38) 



which gives 



By putting Xc = a^fi into Eq. ( 122b and keeping terms up to 0{li), we find 

(39) 

which clearly satisfies the bound Eq. dlSb . and shows that this bound becomes an equality 
when jU — 5- 0. 

The approximation for can be significantly improved by matching the approximation 
for small fi with the behavior of when fi is close to fic- Since rc becomes infinite at 
^ = He, we write 

rc=f /^^,^ , with p = 2^^ (« + VOjS) -— . (40) 



I-M/Mc' ''2/ V / ^x, 

where p is determined such that the correct leading behavior is reproduced when ^ is small. 
The specific form of the divergence at /i = is motivated by the behavior in the case of 
symmetric fitness peaks (see Ea.il47b). 

In Fig. [3] the exact values of rc obtained from numerical calculations are compared with 
the two approximations Eqs. ( 139b and ( 140b . One can find more such graphs in the Mathemat- 
ica file in online supplement. As expected, both expressions give a reliable prediction when 
jU is sufficiently small. However, as the exact value of rc becomes larger, Eq. ( |39b does not 
give an accurate estimation, which is not surprising. The more surprising observation is that 
Eq. ( 140b gives an excellent prediction for rc in almost all ranges. However, Eq. ( 140b becomes 
a poor guidance for rc as t gets smaller, which suggests that the case with small t should be 
treated separately. In the next subsection, we will study the two-locus model for small t. 
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a s = 0.002 and f = 0.04 b s = 0.04 and i = lO"* 




/i fi 

Fig. 3 Comparison of the exact with the approximate solutions Eqs. j39t and 140) for three cases: (a) s 
and (b) t <C s. Eq. 140) shows an excellent agreement with the exact when t is not too small compared to s. 
Eq. )39t is generally poor in predicting r^. 



3.4 Approximation for small f 

Now let us move on to the case with small t {t <^ ^,s) w hich connects the present study to 
the problem with symmetric fitness peaks considered bv lHiggsl ( fl998h . As before, we will 
find Xc from Eq. ( 123b . 

As shown by Hi ggsl lll998l) (see also Appendix |B]|, Xc should approach zero as r — )• 0. 
This can be rigorously shown from Eq. i l23l l. If we assume that Xc is finite as f — s- 0, the 
leading order of Eq. ( 123b becomes 

2s^{2-s){l-2fi) [{(l-2M)x-2(/i,o-/i)}'+4(l-i)/i^o]-«^ = 0, (41) 

where 

^-=2(2^ ^42) 

is the critical mutation probability for the symmetric problem derived in Appendix |B] Since 
the terms in the square brackets are strictly positive, the only real solution is x = 0. 

It might seem natural to assume that Xc = cit + 0{t^) as in Sec. 13. 31 However, this turns 
out to be wrong. Not to be misled by an incorrect intuition, let us expand H{xc) only with 
the assumption that Xc is small, i.e., we do not specify how small Xc is compared to Then, 
to leading order, the equation H{xc) =0 becomes 

2x1 + ^4' - 2a,x,.t^ -a,P = 0, (43) 

where 

a, = (44) 

' 2^(1-2m)(2/|2 + Mco(^-4m))' ^ ' 

Note that terms of order jcj and are neglected compared to x^.. Actually, there is a term 
of order in H{xc), but its coefficient is O(f^), so it is neglected compared to Xct^. Let us 
assume that jc^ ~ . If ^ > 1, the solution of Eq. ( 143b is x^ = —t/2 which does not fall into 
the range < jc^ < xi. If ^ < 1, the leading behavior of Eq. ( 143b should be x^ — a,t^Xc = 
which gives Xc ~ (a/f^)'/^. A more accurate estimate of Xc is derived in the Mathematica 
file in online supplement, which reads 

x,«(a,f2)i/3_ f (45) 
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Note that the power | was aheady observed when we calculated jic for f <C i in Eq. 02b . 
From Eq. i22[ along with Eq. ( 145 l l we get 

J^c « '^cO IH — T7 ^ (flff ) 97 zt , (46) 

where rco is the critical recombination for r = given by (see Appendix IB] for derivation) 

rcQ = , r . (47) 



We will use the same technique as in Sec. l3.3l to improve the quality of the approximation 
for Kc- Since rc diverges at jU = jUc rather than at /icO (note that jUcO > Mc)^ we can associate 
the behavior for small Vc with that for larger Vc in such a way that 

-(l+pof2/^+piO, (48) 



(1-2m)(Mc-M) 



where the coefficients po and pi are determined by requiring that the leading behavior of 
in Eq. l l48t is same as that in Eq. ( 146b . This yields 



Pi) 



where we have used a more accurate expression for /ie than Eq. M2\ . derived in the Mathe- 
matica file in online supplement, which reads 

„ -„ 3(1 -.y) .2^1/3 , 2Mco(1+^) . ™ 

Mc-McO-jT^ r(2McoO H 5 ^- 

For completeness, we present the corresponding which reads 



^"n 4 j -2- '''^ 

Figure|4]compares the exact values with the approximations Eqs. ( 146 l l and ( I48l l. .? and f 
in Fig.|4^ are same as those in Fig.[3l3. For sufficiently small t, both approximations show a 
good agreement. As anticipated, Eq. ( 146b becomes worse as t increases, even though Eq. ( 148b 
is still in good agreement with the exact solution. Needless to say, Eq. ( 148b fails when t is 
not much smaller than s. Although Fig. |4] seems to suggest that the agreement is good even 
for very small pL, this is an artifact because itself is too small. In this regime, one should 
use the approximation developed in Sec. 13.31 

To summarize our findings up to now, we have provided approximate expressions for rc 
which are valid in the specified regimes. Taken together, these expressions cover essentially 
the full range of biologically relevant parameters. 
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a s = 0.04 and ( = 10"^ b s = 0.04 and t = 0.002 




0.005 0.01 0.005 0.01 

Fig. 4 Comparison of the exact with the approximate solutions Eqs. (46) and (48) for s = 0.04 and (a) 
« = 10"'', (b)( = 0.002. For comparison, we also draw the critical recombination probability when t = (rco). 
For t = 10"^, both approximations are in good agreement with the exact solutions. As t increases, Eq. )46t 
starts to deviate from the exact solution, but the improved approximation still has predictive power 



3.5 Frequency distributions 

So far we have investigated the critical recombination and mutation probabiUties. To com- 
plete the analysis, we need to determine the frequency distribution at stationarity. For the 
LFS, this can be readily done using Eq. ( 121b . From Eq. i ll8b we see that the solution with the 
smaller x > will confer the larger mean fitness. Let (xi) be the smaller (larger) positive 
solution. Equation ( 119b shows that B{xs) < B(xi) and A{xs} > A{xi). Note that we are treat- 
ing A and B as functions of x. If we rewrite /3 as 1/ {2^/A/B + 1 +A), it is clear that (xj) 
must be smaller than fi{xi). Introducing 4/3 = h{xi) - f^ix,) and 4/o = Mxi) - fo{xs), 
the mean fitness difference can be written as 

Aw = w{x2) — w{xi) = (wo — wi)A /() + (w3 — wi)zi/3. (52) 

Since Aw is negative and Zi/3 is positive, Afo must be negative. That is, the solution with 
the smaller jc > should confer the larger value of /q. Below we will argue that the solution 
with the larger /o is stable. Hence we limit ourselves to the study of the solution with smaller 
X. We also present approximate expressions for the HFS. 

As before, we have to conduct separate analyses depending on which parameter is the 
smallest. For these analyses we will use the expansions Eqs. ( 139b and ( 146b for rather 
than the improved approximations Eqs. ( 140b and ( 148b . which are not suited for a systematic 
perturbative solution. 

We begin with the case when ^ s,t. The functions h"{xc, r^) and hi (xc) in ( 121b then 
become 

-h"{xc,rc) ^ItP, hi{xc) ^ Oi^s^t^, (53) 



which gives 



.,2 c 



1/2 



S £ 

E, W -Xc I ' I = -XcE, (54) 

(l-f)V«jSM, 

where Xc = a^i^ and the definition of e is clear. If we set x = Xc + £x and r = rc + £r = 
t + Ci^li + Er from ( I39I 1, we get 

t t a 

A=l + -K— -, B=—, ^, (55) 

X Xc{\-E) t[Er-2E_, + {c^-2a^)^i) 
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which are large. Note that A becomes negative when e > 1, which means that the regime 
of validity of the above approximation is quite narrow. To leading order, the population 
frequencies are then obtained from ( 116b as 

\ 1 \ — t FFT 

(56) 

and, by normalization, /o = 1 — /3 — 2/. The second (unstable) solution is obtained by setting 
e -e. 

To find the HFS, we shift the variable x to y = — (x + r) and look for a solution of 
= h(^—t — y) = 0. As shown in the Mathematica file in online supplement, the HFS is 
located at y = (x = —t) for /I — >■ 0. This is a consequence of the fact that when /I = 0, 
the stationary fitness of the HFS is vv = W3. If we now set y = L„>iy«/i" and expand g{y) 
as a power series in \x, the equation g(y) =0 gives (see the Mathematica file in online 
supplement) 



r{\-s-tf ^{s^t){l-s-t) t{s + tf{\-r)^r{ls + t) 



Up to 0{li ), the genotype frequencies for the HFS become 

/«^(l-y2M), /o«A«^/i', (58) 
s + t t 

and/3 = I-2/-/0. 

One can see qualitative differences between Eqs. | |58I > and i l56l l. First, the frequency of 
the less populated fitness peak genotype is different in the two cases. In Eq. ( 156b . /3 = 0{p.), 
but in Eq. ( 158b . /o = 0{)1^). However, this tendency cannot persist when r is large. For 
example, if r = 1, Eq. i llOb suggests that should be 0{p.^) provided / is still 0{)i). Hence, 
this qualitative difference only occurs when r is close to r^. Second, the leading behavior of 
the frequency / of valley genotypes does not depend on r in Eq. ( 158b , which is not true in 
Eq. ( 156b because of the dependence on e. 

To analyze the stability of the solutions, we linearize Eq. ^ at the steady state fre- 
quency. For the stability analysis, we assume that /i(t) = /2(t) for all T, which is true if 
they are equal initially. The linearization then yields a square matrix with rank 2, whose 
largest eigenvalue (in absolute value) determines the stability. For the HFS, the eigenvalues 
up to 0(;ix") are \ — s — t and [l — — r) which are smaller than 1. Hence, the HFS is 
always stable. At r = r^, the largest eigenvalue for the LFS is expected to be 1. Since we are 
restricted to an approximation up to 0{^), all we can show is that the largest eigenvalue of 
the LFS is 1 + O(m^) at r = r^. In the Online Supplement, we show that the largest eigen- 
value for s = t becomes 1 -|- 0{}X^) and the smaller one is { \ — s — t) /{ \ — t). When treated 
numerically, it is easy to see that the stable solution indeed corresponds to the smaller x 
(details not shown). 

The next step we will take is to find the frequency distribution of the LFS in the case 
that t is much smaller than s and p.. From Eqs. ( 145b and ( 146b , we get (up to leading order) 

_^ ^ 6..(2M^ + M.o(.-4M)) . ./3^ 

Mco-M 

hx{x,) = 2s{2-s){\- 2m) (Mco - II ) {a.t^f^ , (59) 
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thus from Eq.( l21b 



-XcT], (60) 



where we have kept Xc = {a,fiyl^ up to leading order from Eq. ( 145b and r] has an obvious 
meaning. Accordingly, A and B become 

A«l+A,, (61) 



where 



The above approximation is valid only when r] ^ 1 (e,. <C f'l^). Note that unlike the previous 
case, A is close to 1 {A, ~ f '/"'). Hence the frequency distribution for the LPS becomes 

f^E(i-2(\-^\ (b,x,+B,£,+Ba + ^\], (64) 



where we have kept the leading order of each term. Since the above approximation requires 
that e,- = r — Fc <^ t^l^ , it cannot reproduce the symmetric solution in Appendix IB] which 
applies when r — ^ at fixed r. 



3.6 Landau theory 

In this subsection, we develop an approximation that is valid when r is close to and the 
asymmetry of the fitness landscape is small, in the sense that t is smaller than all other 
parameters. This approximation is inspired by the Landau theory from the physics of phase 
transitions, and it will allow us to represent both the LPS and the HPS in a simple, compact 
form. 

We start from the observation that, according to Eqs. ( 164b . ( 165b . and ( 166b . the valley 
genotype frequency / w jU in the regime of interest, with the peak frequencies fj, and /o 
symmetrically placed around 1/2 — /i/iw 1/2 — f. Moreover, the difference /o — /a « A; ~ 
r'/^ becomes small for f — >^ 0. This motivates the parametrization 



/0= (^-/)(l-«),/3=Q-/)(l+") 



(67) 
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which defines the new variable u. Inserting this into Eq. (O with = fi we obtain 

w=(l-2/i)(l + (l-M)^). (68) 

On the other hand, from the definition i ll2b of w, we find a relation between / and u such 
that 

, (l-„)(l+,-2^) 2m 
2m 2i + /(l+M) 2s + t{l+u)' 

Up to now, everything is exact. Note that when m ^ 1 and t <^ u, the leading behavior 
of Eq. l |69l l is ^/s which is consistent with the LFS frequency distribution in Eq. l l64b . 
Moreover, as the mean fitness for the HFS in the case of small t is not expected to deviate 
much from 1 — 2jU, the HFS also requires that ? <C m. So for all solutions, the leading behavior 
of / is fi/s. This is rather different from the case when /i is the smallest parameter. 

By keeping the leading terms under the assumption that t <^ ^ <^ s <^ I and t <^ u<^ 1, 
from Eqs. ((Tjl and l |69l l we obtain the equation 

t-iro-r)u-ru^ =0 (70) 

for M, where fq = S/i^/^- If we interpret r as the (inverse) temperature, t as the external 
magnetic field, and u as the magnetization, this has precisely the form of the Landau equation 
for the para- to ferromagnetic phase transition (Plischke and Bergersen, 200©. 
The general solution of Eq. ( 170b can be written in a compact form. Let 

denote the discriminant of Eq. ( l7Qb . When A >0, there is only one real solution which reads 

For r sufficiently far below ro, in the sense that ro — (r^r)^/^, this reduces to 

"HFS w — - — , (73) 
ro-r 

which is the solution of Eq. ( 170b with the cubic term omitted. When A <0, there are three 
real solutions 

MHFS = 2 — — COS u = -2 \ — — sm - T - , (74) 



, 3r 7 3 V 3r ; V6 3, 

where tan0 = 2r^\A\/t with < < n/2. The stable LFS corresponds to the smallest 
value of M, which yields the larger /o among the two solutions with negative u (see the 
discussion in the beginning of Sect. 13.5b , 

«LPS = -2(^)"'sin(f + ^). (75) 

One can easily see that for f — s- (9 — )■ 7r/2) the solutions ( 174b and ( 175b approach the sym- 
metric peak solutions 



mhfs = a/I ~ '"«/'"> "LFS = - a/1 -ro/r. (76) 
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0.001 0.0012 0.0014 0.0016 0.0018 0.002 

r 

Fig. 5 Plots of H = (/s — /o)/(l — 2/) obtained from the exact numerical solutions (symbols) and from 
Eq. )70) (full curves) as a function of r for t = 10^*, s = 10^^, and jl = 10^^. For these parameters ro = 8 x 
lO"-* and « 1.325 x IQ-^. The filled squares are stable solutions and open circles are unstable solutions. 
The dashed line shows the symmetric peak solutions )76t for ryr^o- The approximate solution is seen to be 
valid well beyond the regime where 7] < 1 . 

The critical recombination probability can be found from 4=0 which gives 

- 4 i fe)" ^ 

Note that this agrees with Eq. ( 146 l l only up to order f^/-^ . 

Although the LPS in Eqs. ( 164b . ( 165b . and ( 166b is valid only when Er <S t^^^, the approxi- 
mate solutions Eqs. ( 172b and ( 174b turn out to be in good agreement with the exact solutions, 
provided ro is replaced by the exact critical recombination rate rco for the symmetric peak 
problem [see Ea.(l47b1. In Fig. [5] we compare the exact values of u with the approximate 
solutions for t = 10^*, s = 10^^, and jU = 10^^. For these parameters, r\ becomes larger 
than 1 when w 3 x 10"^. 



3.7 Behavior of the mean fitness 

We are now prepared to discuss the dependence of the mean population fitness w on the 
recombination rate. Since mean fitness is linearly related to the auxiliary variable x through 
dlSb . this amounts to examining how the solutions of ( 120b vary with r. Solving ( 120b for r we 
obtain 

Kx) = -544 (78) 

and r'{x) = H{x) //?, {xf, where H{x) was defined in ( 123b . Since H{x) has a unique root Xc 
in the regime of interest, we conclude that r{x) displays a minimum at x = . Recalling that 
the stable LFS corresponds to the smaller of the two solutions of ( 120b with x > 0, it follows 
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Fig. 6 Behavior of mean fitness as a function of r for (a) ttie HFS in linear scales and (b) tiie LFS in semi- 
logarithmic scales. The parameters are ? = i = 0.1 and = 0.001. For the purpose of clarity, we show the 
difference between w and its limiting value for large r, which is W3(l — 2^11) for the HFS and wq(1 — 2jl) 
for the LFS. (a) Mean fitness decreases with r. (b) Mean fitness increases (decrease) with r for the stable 
(unstable) LFS. Note that this figure displays the positive quantity wq — m>/(1 — 2/i). 



that the mean fitness of the stable (unstable) LFS increases (decreases) with r. In addition, 
the fitness of the HFS must be a monotonic function of r. The behavior of all three solutions 
is illustrated in Fig|6l which shows that fitness decreases with r for the HFS. 

These results can also be deduced from the approximate solutions given in Sects. 13.51 
and l3.6l In particular, taking the derivative of ( 168b with respect to u we see that dw/du < 
always. Since ii is an increasing (decreasing) function of r for the HFS (LFS), it follows that 
fitness decreases with r in the former case but increases in the latter. For ?- — >■ oo the solutions 
| |74I75| | approach mhfs — ^ 1 and mlfs — ^ — L respectively, with corresponding limiting fitness 
values vi>HFS ^ >V3(I — 2fi) and wlfs h'o(1 — 2jU). As can be anticipated from Eq. ([8]l 
(see also discussion in Sect. 13. lb . the limit is approached from above for the HFS but from 
below for the LFS. Note that in the symmetric case (t = 0) the fitness is w = m'3(1 — 2/x) = 
wo(l — 2jU) independent of r for r > (see Appendix IbTi. 



3.8 Asymmetric valley fitnesses 

In this final subsection, we will consider briefly how the results would be affected if wi = 
I —t — si and W2 = 1 — f — J2 with si =^ S2 (without loss of generality, we can set si < S2). 
We first show that bistability requires reciprocal sign epistasis, i.e. both and S2 have to 
be positive. To see this, suppose that Ji < and ^2 > 0, such that the ordering of the fitness 
values is W2 < wo < wi < W3. Then positivity of Eq. ^ requires w > (I — 2jU)wi or vi> < 
(1 — 2jU)w2 < W2. The latter possibility is ruled out because the mean fitness cannot be lower 
than the fitness of the least fit genotype, and the former inequality contradicts the condition 
w < ( I — 2fi)wQ imposed on the LFS by the positivity of Eq. We conclude that only the 
HFS with mean fitness in the range (1 — 2/i)w3 < w < W3 can exist. 

To extract some information about the case S2 > si > 0, we introduce the variable m in a 
similar way to Eq. ( 167b such that 



/0 = (1-/1-/2)^, /3 = (1-/1-/2)^ 



(79) 
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which yields again Eq. i l68t for w. The above parametrization along with Eq. ([9]l gives 

2u{s2-si) 

J\=J2 + 777- — — h- (oO) 

Thus, if .52 — ii ^ i2> setting /i = /2 = / is not a bad approximation and the results 
presented above remain valid. 

When S2 — si is comparable to the other parameters the calculation is much more com- 
plicated. We will not treat this case in any detail, but we can provide necessary conditions 
for bistability. First, the necessary condition for r given in Eq. ( 115b is still valid, because 
Eq. ( 115b is obtained only from the denominator of Eq. dlOb . A necessary condition on /i 
similar to Eq. ( 114b is also available from the requirement that W2<w < (1 — 2jU)>vo, but the 
result is not very useful because it is independent of the magnitude (or sign) of ii. To get a 
refined condition, we need some futher analysis. 

From the definition of vt>, we find (see also the Mathematica file in online supplement) 

. {t{2li{\-u)+u^-\)+Aiiu){u{s,+S2)+t{u+\)) 

u{si{4s2U + t{u+l)2)+t{u+l)^S2+t)) ' ^ ' 

which should be smaller than 1 for —1 < u < 1. Now assume that a critical value exists 
(rc < r < 1) beyond which bistability appears. For r = I Eq. dlOb shows that /0/3 = /1/2 
for any equilibrium solution. On the other hand, we expect that for small ^ the frequencies 
of both valley genotypes are small, /1/2 ^ 1. Computing /o/? using Eq. ( 179b this is seen to 
imply that u for the LFS should be very close to — 1 when r= I. Expanding Eq. JSlb around 
M = — 1 gives 

/j+y2=^(£l±£iKi^ + ^(,j+,2-M(2 + >?i+^2-2f))(l+M)+0(l+M)2. (82) 
For a LFS to be possible, the leading order term should be smaller than unity, which gives 

^ (^i+^2)(l-0" 2(1-0' 

where sh is the harmonic mean of si and ^2. Note that Eq. ( 183b reduces to Eq. ( 114b when 
si =i2-Hence, ifO< Ji ^J2>Mrnustbe smaller than /(I —f) to have multiple solutions. 



4 Discussion 

In this work we have presented a detailed analysis of a deterministic, haploid two-locus 
model with a fitness landscape displaying reciprocal sign epistasis. We have established 
the conditions for the occurrence of bistability, and derived accurate approximations for 
the critical recombination rate at which bistability sets in. For r < there is a single 
equilibrium solution in which the fittest genotype is most populated. For r > r,; we find two 
stable equilibrium solutions, one of which is concentrated on the fittest genotype (the HFS) 
and a second one concentrated on the lower fitness peak (the LFS). For jU — )• these two 
solutions become point measures, in the sense that /; — s- 1 and /o — )■ 1, respectively, but for 
any finite mutation probability all genotypes are present at nonzero frequency. 

We briefly summarize the most imporant quantitative results presented in this paper. The 
expressions ( |39b and ( 146b for are based on a systematic expansion in terms of the muta- 
tion probability jU and the peak asymmetry t, respectively, while the interpolation formulae 



20 



Eq. i l40b and Eq. ( 148b provide numerically accurate values of r^. over a wide range of param- 
eters. In particular, our results show that the lower bound ( 115b on Vc becomes an e quality for 
/i ^ 0, which is consiste nt with earlier results obtained either directly at /X = iFeldman , 
197 1 ; Rutschman', 'l 994h or using a unidirectional mutation scheme dEshel and Feldmaii , 
~970. : .Kariin and McGregojJ 197 Ih . Clearly the limiting behavior for jU ^ should not de- 



pend on the mutation scheme employed. Approximate results for the stationary frequency 
distributions are found in Eqs. ( 156b . i l58b and Eqs. ( 164b . i l65b . i l66b . Of particular interest are 
the simple formulae derived from the cubic equation i l70b . which are remarkably accurate 
close to the bistability threshold and for small t. 

The consequences of our results for the question of a possible evolutionary advantage of 
recombination are twofold. Dynamically, the onset of bistability implies that recombination 
strongly suppresses the escape of large populations from suboptimal fitness peaks. In the 
deterministic infinite population limit considered here, the escape time diverges at r = 
l^ain, 2010), whereas in finite populations one expects a rapid (exponential) incr ease of the 
escape time for r > dStephanl 1 19961 : iHiggsi 1 19981 : IWeinreich and Chad I2OO5I) . In a mul- 
tipeaked fitness landscape, recombination can therefore dramatically slow down adaptation 
lldeVissereta]L l20091. On the other hand, we have seen in Sect. 13.71 that the equilibrium 
mean fitness may increase or decrease with the recombination rate when r > depending 
on which of the two equilibria is considered. In a fitness landscape with more than two 
peaks one anticipates an even richer structure of stationary solutions with a correspondingly 
complex dependence on recombination rate. 

It would be of considerable interest to extend the present study to fi nite populatio ns. For 
the case of symmetric peaks {t = 0) this problem has been addressed bv lHiggsl ( fl99i) in the 
framework of a diffusion approximation. A key step in his analysis was the reduction to a 
one-dimensional problem by fixing the frequency of the valley genotypes at its stationary 
value f = p. /s. However, we have seen above that in the case when s and t are large in 
comparison to /I, / cannot be treated as a constant and therefore the reduction to a one- 
dimensional diffusion equation is generally not possible. Some progress could be made in 
the regime where the (effectively one-dimensional) Landau equation l l70b applies, and we 
intend to pursue this approach in the future. 
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A The function h{x) and its values at selected points 

All information in tliis appendix can be found in the Mathematica file in online supplement, and is provided 
here only for completeness. 

By squaring both sides of Eq. 1171 . the steady state mean fitness becomes the solution of the equation 
h{x) = where 

.1 B , , , A{i-t-x)h(x) 
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As shown in the Mathematica file in online supplement, h{x) takes the form ho{x) + rh[ (x) with 

ho (x) = b4X* + b3X^ + b2X^ + bix-ho, (85) 
hi{x) = -(l-2ti)^C3X^-{l-2ii)c2X^+cix-co, (86) 

where 

b4 = (1-2a()(2s + «), C3 = 2.s + f-(^ + f)^ 

b3 = {t^ + 2st-2s^){l-2ii) + n^{4c3+t^), Q = (/ + 2^1 - 4<;i)c3 - 
b2 = -3s^t{l-2ll)-ll^ [4{l-2t)c3+3t^{l-t)] , 
ci = f(l-2^l)(/-2Al(l-0c3) + AlV(l-.v-f)2, 
b, = -(l-2Ai).vV-/J^r[(4-5r)c3-?(l-?)(2-30], 

CO = (l-t){l-s-tft^^^, bo = {l-t){2-s-2t)st^^^. (87) 
Note that b4> if fl < ^. The values of /i(jc) al x = 0, x = —t, and A' = jri = wq — / {i —2n) are 

/i(0) = -foo-rco, (88) 
H-t) = -fV[l-(l-'-)(l-.v-f)'], (89) 
M-vi) = - (T^lljT [1 - (1 - '■) (1 - + - Ai) - (90) 

1 

2- 



B Solution for symmetric fitness peal^s (f = 0) 

When f = 0, our problem is reduced to that (approximately) solved bv lHiggj il998l) . For the paper to be 
self-contained, we solve the case with t = exactly in this appendix. 

When 1 = 0, Eq. {8) suggests either = /o or vv = 1 — 2^11. Let us first consider w = 1 — 2/1. Needless to 
say, this solution is impossible if /I > j . Since v? = (/o + ) + 2/( 1 — i) = 1 — 2fs, we get 

/=-,/o + /3 = l-2^. (91) 

s s 

From Eq. jlOt with w = 1 — 2jl, we get 

(/0-/3)' = (/o + /3)'-4/o/3 = ^(l-2/i)^, (92) 

where ^ = (2 — s) (flco — /J ) (r — r^o) with r^o and fl^o given in Egs. j42t and (47). 

To have an asymmetric solution (/o ^ /3), ^ should be nonnegative. Because /I > Hco implies r^o < 0, 
^ is nonnegative only if r > r^o and n < Hco (note that when n = flco, ^ is nonzero and negative, because 
Eq. 1471 diverges as 1/ {fl^o — /J) for fl — > /ico)- Since r^o cannot be larger than 1, the more restrictive condition 
on n becomes li < jXc where jic is the solution of the equation r^a = 1 given by 

M. = ^. (93) 

Note that /i^ is the same as /icO up to leading order of .v. 

Now let us find the solution with f3 = fo- From 1 — 2fs = w and Eq. jlOt along with the substitution 
w = (1 — 2/i)(l +y), we obtain 

g,(>.)^/ + (r(l-,v)+. + ^))' + ^ =0. (94) 

Since 



g,,(-^) = -(l-.)(2-,v) 



-r{ii + llco) 



1-2/i 

is negative, there is only one solution with y > —s which is 

y= — nr- (95) 

When 1^ < (either jl > fl^o or /I < Hco together with r < r^o), the larger solution is nonnegative. On the 
other hand, if ^ > (p < Hco and r > rco), the larger solution which is still larger than — .s is negative. 
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